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Abstract 

The  invariant  torus  of  a coupled  system  of  Van  der  Pol  oscillators  is  approximated  using 
bicubic  B-splines.  The  paper  considers  the  case  of  strong  nonlinear  coupling.  In  particular,  the 
shapes  of  invariant  torii  for  the  Van  der  Pol  coupling  parameter  A are  computed  in  the  range  [0.1, 
2.0].  Comparisons  are  given  with  results  obtained  by  the  MATLAB  differential  equation  solver 
ode45.  Very  good  normed  residual  errors  of  the  determining  equations  for  the  approximate  tori 
for  the  cases  A = 0.1,  0.6  are  shown.  At  the  upper  limit  of  A = 2.0  memory  errors  occured 
during  the  optimization  phase  for  solving  the  determining  equations  so  that  full  optimization 
for  some  knot  sets  was  not  achieved,  but  a visual  comparison  of  the  resulting  invariant  torus 
figure  showed  a close  similarity  to  the  solution  using  ode45. 

Keywords:  bicubic  B-splines;  determining  equations;  invariant  torus;  large  parameter  case; 
optimization;  Van  der  Pol  oscillators 


1 Introduction 

Physical  problems  involving  coupled  vibrating  systems  often  lead  to  ordinary  differential  equations 
of  the  form 


x + Ax  = AN(x)  (1) 

where  x 6 Rn,  A is  diagonal,  and  N is  a nonlinear  function  of  x.  A is  a parameter  that  scales 
the  nonlinearity.  There  are  many  forms  of  solutions  to  systems  like  (1),  but  we  are  interested  in 
approximating  a particular  class  of  solutions.  In  particular,  we  are  interested  in  systems  that  exhibit 
an  invariant  surface  of  solutions.  That  is  a surface  on  which,  if  a solution  passes  through  a point 
on  the  surface,  then  the  entire  solution  remains  on  the  surface.  In  the  case  of  periodic  solutions 
the  surface  is  sometimes  referred  to  as  an  invariant  torus.  In  the  case  of  small  parameters  there 
is  a large  literature  on  the  existence  and  approximation  of  invariant  tori  (see  e.g.  Bogoliubov  and 
Mitropolsky  [1],  Hale  [9.  10],  and  Gilsinn  [5]).  For  the  case  of  a large  parameter  there  is  not  much 
in  the  way  of  literature.  The  work  of  Ge  and  Leung  [4]  showed  that  the  Fast  Fourier  Transform 
could  be  used  as  a tool  to  approximate  invariant  tori  in  the  large  nonlinearity  case.  In  this  study  we 
introduce  a bicubic  B-spline  computational  algorithm  that,  at  least  for  a certain  class  of  nonlinear 
systems,  computes  invariant  tori  to  a high  degree  of  accuracy  for  a larger  range  of  parameters.  In 
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particular,  following  the  work  of  Gilsinn  [6,  7],  we  shall  consider  the  system  of  coupled  Van  der  Pol 
oscillators 


x\  + iqx  i = A(1  — x2  — cv  i x 2 ) x i , 

X~2  + ^2x2  = A(1  — 02%!  — x\)x,2,  (2) 

where  ^1/^2  are  irrational  multiples. 

Under  the  change  of  coordinates  x\  = r\^2  sin(/xia;i),  x\  = //iri1//2  cos(/.iiu;i),  x2  = sm(fi2^>2), 

and  £2  = LLor-21/2  cos(/i2Ui2)  this  system  becomes 


cJi  = 1 — (—  A f ) (sin(/.iiu;i)  008(^10;!)  — ri  sin3 (^iicui)  cos(/.iicui) 

V hi  ) 

—air-2  sin(/.iicui)  cos(/iiu;i)  sin2 (^2012)), 

cJ2  = 1 + A ( -—  ) (sin(/U2u;2)  cos(/i2c<;2) 

V h2/ 

—0:2X1  sin2 (/xiuii)  sin(/X2^2)  cos(^2^2)  - r2  sin3(/i2cu2))  cos(p.2ui2),  (3) 

r'i  = 2Ar!(cos2(/riu;i)  — ri  sin2(/uia;1)  cos2(/riu;i) 

— Oi7'2  COS2(/7itUi)  Sin2(/U2CU2)), 

r'2  = 2Ar2(cOS2(/i2^2)  — 02Xi  sin2(//itUi)  COS2(/72^2) 

— £2  sin2(//2<£'2)  cos2 (1x2^2)), 

which  is  of  the  form 


S2  = d + A©(f2.  r), 

f = AR(0,r), 


with  f2  = 


UJl 

U>2 


(4) 


©(f2,  r) 


/ ©i(fi,r)  \ 
V ©2(fi,r)  / 


^ (-^-)  (sin(/r1tu1)cos(/.iiu;1)  ^ 

— ri  sin'!(/uiu;i)  cos(/itiu;i) 

-Oi\r2  sin(/iiu;i)  cos(/r1u;1)  sin2 (^2^2)) 

(_i)  (sin(h2W2)  cos(/x2^2) 

-o2ri  sin2(/iitui)  sin(/,i2^2)  cos(/r2ai2) 

\ -r2  sin3(/r2W2))  cos(/x2^2)  / 


R(T2,  r) 


V R2(^,r)  yi 

/ 2ri(cos2(/u1u>i) 

— ri  sh\2(niuJi)  cos2 (hiuji) 
-Oir2  cos2(/ti1cui)  sin2 (^2^2 )) 
2r2(cos2(//2w2) 

—02X1  sin2(^io>i)  cos2 (+12UJ2) 
„ „ 2 / . , . . \ 2 1 . . . \\ 


\ 


/ 
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For  the  current  study  we  will  assume  that  an  invariant  torus  exists  for  (3)  in  the  form 


S = {(/xicJi,/i2W2,n,r2)  : n = n2U2),r2  = H2U2), 

(Fi^i,  H2U2)  € T~ } . (5) 

Using  the  invariance  of  the  graphs  of  /(piuq,  ^2^2)  and  g(niuji,  H-2^2)  under  the  dynamics  gen- 
erated by  (3)  we  have 


r = DLJiF(hilj1,H2U)2)ui  + 

= AR(n,F(/iiWi,/i2w2)), 


where  F = 


/ 

9 


Now,  from  (4), 


(6) 


oil  — 1 + AOi  (fi,  F(pic<;i,  p2<^2))) 

u>-2  = 1 + A02(fi,  F(^icai,  /i2cu2))-  (7) 

We  therefore  substitute  (7)  into  (6)  to  get  the  following  quasilinear  partial  differential  equation  that 
F^iuq,  H2U2)  must  satisfy  in  order  for  its  graph  to  be  an  invariant  torus  (See  [8],  [13]). 


N {¥{ni'jjx,  112^2))  = 
+ 


For  the  remainder  of  the  discussion, 
notation  in  the  equations. 


£UiF(piu;i,  ^2^2) (1  + A©i  (fh  F(pitui,  H2U2))) 

F(pitui,  p2^2)(l  + A02(f2,  F(jiiWi,  /i2cj2))) 

AR(0,  F(/iicui , /r2tJ2))  = 0.  (8) 

we  make  the  following  substitutions  in  order  to  simplify  the 


6 — Hi  ui\, 

<P  = H2^2- 


2 Bicubic  Spline  Collocation  Method 

We  now  use  the  method  of  collocation  with  B-splines  to  approximate  the  functions  rq  = f{9 , <p)  and 
r2  = g{0,  (p).  We  let 

0 = 90  < 9i  < . . . < 9n  = 2t r, 

0 — 00  (pi  . . . <C  (pm  — ^71, 

be  a uniform  partition  of  the  rectangle  ranging  over  [0,  27t]  in  the  9 direction  and  [0,  2tv]  in  the 
(p  direction  where  n and  m are  the  number  of  knots  in  the  9 and  <p  directions,  respectively  and 
A h = 9j+i  — 9j  = 27 r/n,  i — 0, . . . , n — 1 and  A k = (pj+i  — (pj  = 2ir/m:  j = 0, . . . , m — 1.  We  seek 
bicubic  spline  surfaces  of  the  form 

71  + 1 771+1 

ri,nm(e,<p)  = J2  *jBi(0)BM, 
i=— 1 j=—  1 

7i+ 1 ra-fl 

- E E dlJB,(9)Bj((P),  (9) 

i=— 1 j=-i 
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where  {B-i(9),  Bq(0),  . . . , Bn+i(0)}  are  the  B spline  basis  functions  in  the  0 direction  with  knots 
Oh,  h = 0, . . . , n and  {5_i(0),  £?o(0), . . . , Bm+i ((f))}  are  the  B spline  basis  functions  in  the  0 direction 
with  knots  <pk,  k = 0, . . . , m.  The  functions  B,(0)  are  defined  in  Prenter[12]  as  follows  by  introducing 
four  additional  knots  in  each  direction  satisfying 

0-2  < 0-i  < 0o  = 0 and  0n+2  > 0n+1  > 0n  = 2n, 

0- 2 < 0-1  < 00  = o and  0n+2  > 0n+i  > 0n  = 27t,  (10) 


*<*>  = *< 


(0-0,- 2)3  if  6>  e [0i-2,0i-i) 

h3  + 3 h2(0  - 0i-x)  + 3 h(0  - 0 1—\ )2  - 3 (0  - 0,-,)3, 

if  0 E [0i- 1, 0,\ 

h3  + 3 h2(0l+1  -0)  + 3 h(0i+1  - 0)2  - 3 (0i+1  - 0)3, 

if  0 G [0i,  0I+ 1] 

0 otherwise. 


(11) 


The  function  Bj(<f>)  is  defined  similarly.  These  functions  are  twice  continuously  differentiable  on 
the  real  line.  We  can  therefore  define  B((0)  and  B"(0)  using  the  equations  for  B,(0).  To  compute 
the  approximations  r\^nrn  and  r2,nm  in  (9)  at  the  knots,  we  note  that  by  substituting  the  knots  Oh 
and  0^  into  the  definitions  of  B , B' , B", 


Bi(0h)Bj(<t>k) 


16  if  h = i,k  = j, 

4 if  h = i — 1,  k = j or  h = i,  k = j — 1 

or  h = i,  k = j + 1 or  h — i + 1 ,k  = j , 

< 1 if  h = i — l,k  = j — 1 

or  h = i — 1 , k = j + 1 

or  h = i + 1 , k = j — 1 

or  h = i I- 1,  fc  — j ~1  1, 


(12) 


and  Bi(0h)Bj(4>k ) = 0 for  h > i + 2,  k > j + 2 and  h < i — 2,k  < j — 2. 

The  partial  derivatives  of  the  interpolants  ri<nm(0,  0),  r2^nm(0,  0)  at  the  knots  can  be  easily  found 
using  the  first  derivatives  and  the  products  B((0h)Bj((f)k ) and  Bi(Qh)B'j(<f)k)  of  the  spline  functions 
in  the  tables  below. 
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(^-ii0j-i)  (0j— i ,4>j)  {6j-i,4>j+\) 


BtmM 

3 

Ah 

12 

Ah 

3 

Ah 

BMB'jW 

3 

Ak 

0 

3 

Ak 

(gf,0j-l)  (flj,0j)  (6>,,0i  + i) 


0 

0 

0 

Bi(0)B'(ct>) 

12 

Ak 

0 

12 

Ak 

($i+l ) *Pj  — 1 ) 

{Oi+i , <pj) 

(#i  + l 5 0j  + l) 

B'i(0)Bj 

(<t>) 

3 

Ah 

12 

Ah 

3 

Ah 

Bi(0)Bj 

(0) 

3 

Ak 

0 

3 

Ak 

Table  1:  Partial  derivatives  of  the  interpolants  riiTl7n(0,  <t>),r2,nm(0,  0)  at.  knots 

Using  the  above  information,  the  values  of  the  approximations  ri,„m  and  r2,nm  at  the  knots  can  be 
expressed  in  terms  of  the  parameters  ctj  and  d,j  as 

U ,nm{fihi  0fc)  — 16C/j + 4 (Cfi—l,k  T Ch,k  — 1 T Qi,fc+1  T Cf i+l,A:) 

+C/i-l,fc-l  + + Ch+pfc-l  + Cft+i^+i, 

^2,nm  (0/u0fc)  — 16d/i,A;  + 4(^-1^  + dh,k—l  + ^/i,fc+l  + dh+l,k) 

4-dh-i,k-i  + dh-i,k+i  + dh+i,k-i  + dh+\,k+i-  (13) 


A similar  computation  shows  that  the  derivatives  m {Oh,  <j>k),  '"'ai'""  {®h,  0fc),  "00’"  (@h,  4>k)  and 


dr  1 


dr  2, 


dr2,n 


jj~(0h , <Pk)  at  the  knots  have  the  following  representation  in  terms  of  the  parameters. 


dr  i,nr 

00 

dritnr 

dcj) 


:(0h,4>k)  = 

■( Oh,<j>k ) = 


Ah 


(Cfc+1,A+1  + ^Ch+l,k  + Ch+  l,fc-l 


— Ch-l,k+\  - ^Ch-\,k  — Ch-l,k-\), 

3 

-^{ch+l,k+ 1 — Qi+l,fc-l  + 4Ch,k+ 1 
— 4c/l)fc_i  — Ch-i,k-i), 
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dv-2  ,nm 

do 


{QhAk) 


df'  2,nm 

<90 


{Oh,  0fc) 


~^j^(dh+i,k+i  + ^dh+i,k  + dh-\-\,k-i 
—dh-i,k+i  ~ ^dh-i^  — dh-i,k-i), 


3 

AA? 


(dh+l,k  + l 


dh+i,k-i  + 4rf/,,fc+i 


—4 dh,k-i  + dh-\,k+i  ~ dh-\,k-i)‘ 


(14) 


To  apply  the  collocation  method,  we  choose  the  collocation  points  to  coincide  with  the  knots  and  then 
substitute  the  nodal  values  ri.nm(#/j,  0fc),r2,nm(#/, 4>k)  and  their  derivatives  obtained  in  equations 
(13)  and  (14)  into  the  partial  differential  equation  (8).  This  results  in  the  following  system  of 
(n  + 1)  x (m  + 1)  coupled  nonlinear  equations  in  (n  + 3)  x (m  + 3)  parameters  Cij  and  di-y 


-^j-(Ch+l,k  + l + ^Ch+iy  + Ch+  i,fc-i  — Ch-l,k+l 


— ^Ch-\,k  ~ Ch- 1 , fc  — 1 ) 


((1  + A@i  (f2,  ri^nm(0h , 0/c)))) 


+ 


^r(cM-i,fc+i  _ Ch+i,k-i  + 4 Ch,k+i  ~ 4 Ch,k-i 


+Ch-l,k+l  ~ Ch-l,k-l) 


(1  + A@2(fi,  ritnrn(6h,  <t>k))) 


-XR(n,rliTim(6h,<f>k))  = 0, 


(15) 


for  h = 0, . . . , n and  k = 0, . . . , m. 

3 Periodicity  Conditions 

The  surface  we  are  approximating  is  a closed  one  generated  by  two  periodic  closed  curves.  The  spline 
model  is  therefore  chosen  so  that  it  is  closed  in  both  the  6 and  0 directions.  This  is  accomplished 
by  repeating  both  the  knots  and  the  control  vertices  as  follows. 

(0n,0fc)  = (0o,<Pk),  0 <k<m,  {Oh  Am)  = (Qh,<j>o),  0 <h<n 


C-l.j  — Cn—ij 

Cn+ij  = Cij,  -l<j<m+l 


Ci,  — 1 — Q,m— 1 

Ci,m+ 1 = C/,1,  -1  < i < n + 1 


Cn,k  = c0,it,  0 < k < m 

Ch,m  = Ch, 0.  o < h < n 

Application  of  the  above  periodic  boundary  conditions  reduces  equation  (13)  to  a nm  x nm  block- 
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pentadiagonal  matrix  equation  of  the  form 


^l,nm(0)  0) 

' D 

B 

0 

0 

B ' 

Co,0 

T l,nm(0,  1) 

B 

D 

B 

0 

0 

Co,l 

ri,nm(  0,2) 

= 

0 

B 

D 

B 

0 

Co, 2 

ri,nm((n  - 1),  (m  - 

2)) 

0 

0 

B 

D 

B 

C(n  — l),(m— 2) 

ri,nm((n  - 1),  (m  - 

!))  . 

B 

0 

0 

B 

D 

C(n— l),(m— 1) 

(16) 


where  D and  B are  m x m pentadiagonal  matrices  with  entries 


and 


16  4 0 ...  0 4 

4 16  4 0 . . . 0 

0 4 16  4 ...  0 


0 . ..  0 4 16  4 

4 0 ...  0 4 16 


B 


41  0 ...  0 1 

14  1 0 ...  0 

01  4 1 ...  0 


0 ...  0 1 41 

_ 1 0 ...  0 1 4 _ 

with  a similar  matrix  equation  obtained  for  the  approximation  r2,nm  by  replacing  the  control  points 
Ci,j  with  dtj. 

Similarly  the  system  of  partial  derivatives  of  ri>nm  in  (14)  reduce,  after  applying  the  periodic 
boundary  conditions,  to  a block  diagonal  form 


%^(0,0) 

r o 

B 

0 

0 

—B  1 

Co.O 

%P(0,1) 

-B 

0 

B 

0 

0 

Co.l 

sT(0, 2) 

= 

0 

-B 

0 

B 

0 

Co, 2 

aem^n  l)^m 

2)) 

0 

0 

-B 

0 

B 

C(n—  l),(m— 2) 

L aP  dn  !)>(m 

1))  J 

B 

0 

0 

-B 

0 

C(n— l),(m— 1) 

(17) 


where  B is  a m x m diagonal  matrix  with  entries 


4 1 
1 4 


0 ...  0 1 

1 0 ...  0 

4 1 ...  0 


0 ...  0 1 4 1 

1 0 ...  0 14 
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and 


1 , nm 

d<t> 

drltnm 

d(f> 


^r(o,o) 

^r(0,i) 

r d 

B 

0 

0 

B 1 

Co.o 

B 

D 

B 

0 

0 

Co.X 

^r(o,2) 

= 

0 

B 

D 

B 

0 

Co, 2 

((n-  1),  (m  — 
((n  - 1),  (m  - 

2)) 

1)) 

0 

B 

0 

0 

B 

0 

D 

B 

B 

D 

1 

? 3 

1 1 

? ? 

1 1 

H-*  to 

1 

where  D and  B are  m x m diagonal  matrices  with  entries 


and 


0 1 0 ...  0 -1 
-10  1 0 ...  0 
0-10  1 ...  0 


0 ...  0 -1  0 1 

1 0 ...  0 -1  0 


0 1 0 ...  0 -1 
-10  1 0 ...  0 
0-10  1 ...  0 


0 ...  0 -1  0 1 

1 0 ...  0 -1  0 


(18) 


with  similar  expressions  obtained  for  the  partial  derivative  approximations  '>'2QQTn , by  replac- 

ing the  control  points  Cjj  with  dij. 


4 MATLAB  solution 

We  now  describe  the  algorithmic  details  involved  in  computing  the  bicubic  B-spline  coefficients 
Cij  and  dij  by  setting  up  an  appropriate  objective  function  to  be  mimimized  using  MATLAB’s 
nonlinear  minimization  subroutine.  Our  objective  function  simply  computes  the  “residual  error”: 
the  difference  between  the  derivative  of  the  approximations  r i'nrn,  r2/nm  and  the  derivatives  f i,  r2 
respectively  as  given  by  the  system  (3)  which  is  equivalent  to  the  partial  differential  equation  given 
by  (8).  Substituting  equations  (16), (17),  and  (18)  into  the  partial  differential  equation  system  (15), 
we  reduce  it  to  a 2nm  system  of  coupled  nonlinear  equations  in  2nm  parameters  Cj  and  d,3.  The 
determining  equations  can  now  be  written  as 

Ri,hk(A)  = m 1 + A0i(n,  r(h,  k))))  + fc)(l  + A02(U.  r (h,  k))) 

-ARi(fi,r  (M))  =0, 

r\  A 

R2m(A)  = k)((l  + A 0X(U,  r (h,  k))))  + k)(  1 + A02(U,  r (h,  k))) 

06  00 
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— AR2(fi.  r (h,  k ))  = 0, 

(19) 

for  h = 0, 1, . . . , n,  k = 0, 1, . . . , m,  where  A = [c;  djwith  c = (c_q_i,  c_ qo,  • • • , Cn+i}m+i)T  and 
d = d-i,o,  • • • , hn+iiTn+i)r.  Equation  (19)  can  then  be  used  to  find  an  approximation  for 

the  vector  A by  a minimization  process,  such  as 

2 nm 

min^|f?fc(A)|2.  (20) 

c.d  — 
k= 1 

We  use  MATLAB’s  f solve  command  to  solve  the  system  of  nonlinear  equations  for  the  parameters. 
For  an  initial  guess,  we  fix  the  following  parameters  for  this  case  study: 


Qi  = 0.2,  q2  = 0.4,  pi  = 1,  p-2  = 1.414,  (21) 

and  take  the  first  approximation  to  the  radii  as 

ri,nm  = 3.53,  r2,nm  = 1-18.  (22) 

We  note  that  the  matrix  form  of  equation  (16)  is 

Ac  = r,  (23) 

with  the  matrix  A being  nonsingular.  That  is,  system  (16)  has  a unique  solution  and  consequently, 

the  collocation  approximate  (9)satisfying  the  periodic  boundary  conditions  is  uniquely  defined.  We 
can  therefore  solve  (23)  by  substituting  the  initial  guess  for  r i>nm  and  r2,nTn  as  3.53  and  1.18  respec- 
tively into  the  right  hand  side  of  (23)to  obtain  the  vector  c of  the  unknown  parameters  ctJ  as 

c = A~lr.  (24) 

A similar  calculation  yields  the  initial  vector  d of  the  unknown  parameters  dij. 


5 Optimization  Issues 

For  the  first  few  runs,  we  coded  the  objective  function  and  called  up  the  f solve  optimization  routine 
without  supplying  the  sparsity  pattern  of  the  Jacobian.  This  worked  well  for  small  values  of  A but 
as  the  strength  of  coupling  increased,  we  had  to  increase  the  size  of  the  knot  domain  for  getting  a 
good  approximation  as  well  as  for  the  optimization  to  terminate  successfully.  Increasing  the  size  of 
the  knot  domain  for  bigger  values  of  A slowed  down  the  optimization  considerably.  Computationally, 
the  model  proved  expensive  even  for  A = 0.5,  which  appeared  to  require  a 40  by  40  knot  domain  for 
a reasonably  good  approximation. 

The  residual  error  for  any  point  in  the  knot  domain  depends  only  on  the  nine  surrounding 
parameters  and  therefore  has  very  few  non-zero  partial  derivatives.  The  Jacobian  is  sparse  with 
well-structured  blocks  of  non-zero  elements.  The  computational  cost  of  finite  differencing,  which 
MATLAB  uses  to  approximate  the  Jacobian,  can  therefore  be  reduced  significantly,  if  we  can  supply 
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the  structure  of  the  Jacobian  to  MATLAB.  This  is  done  by  setting  JacobPatternto  ‘ on’  in  options. 
This  was  our  first  improvement  to  the  code.  With  this  enhancement  the  results  were  encouraging 
for  A between  0.1  and  1.0:  The  number  of  function  evaluations  per  iteration  dropped  significantly 
(for  example,  for  A = 0.5  with  a 40  by  40  knot  domain,  the  number  of  evaluations  dropped  from 
3202  to  just  32  per  iteration!), and  the  time  required  for  the  optimization  to  successfully  terminate 
reduced  considerably.  The  following  table  summarizes  these  results.  Unfortunately,  supplying  the 


A = 0.5,  knots=40 

Implementation 

Normalized  Residual  Error 

Time  for  Optimization 

Jacobian  Structure  Unspecified 

2.26395e-007 

4.0188e+003  secs 

Jacobian  Structure  Specified 

2.33638e-007 

41.3750  secs 

Jacobian  Sparsity  Pattern  appears  to  be  insufficient  for  A bigger  than  or  equal  to  1.0,  which  require 
bigger  knot  domains  for  achieving  the  required  accuracy.  For  these  coupling  strengths,  running  the 
code  generates  an  ‘Out  of  Memory’  output  in  MATLAB.  The  figures  below  for  A = 1.0  & 2.0  were 
generated  with  a 50  by  50  knot  domain  arid  the  optimization  was  incomplete.  Although  unoptimized, 
the  bicubic  spline  surface  appears  to  capture  the  surface  generated  by  numerical  integration  of  the 
full  system  as  shown  in  figures  (3)  and  (4)below. 

6 Results 

We  present  here  the  results  of  running  our  code  for  four  different  values  of  the  coupling  parameter.  In 
each  case,  the  error  trends  were  computed  between  the  numerically  integrated  angle-radial  equations, 
equation  (4)using  the  MATL  AB  ode45  solver,  and  the  approximate  system  of  phase  equations  on 
the  torus,  given  by 

= d + A@(fi,  rnm($,  </>)), 

r = rnm(<9, 0).  (25) 

We  set  the  error  tolerance  in  these  computations  to  10-8  and  tabulated  the  results.  First,  we 
consider  the  case  A = 0.1.  We  ran  the  code  for  three  different  sizes  of  the  knot  domain  starting 
with  20  knots  in  each  of  the  0 and  0 directions.  Figure  (l)shows  the  full  and  approximate  system 
integration  for  the  case  A = 0.1. 


Figure  1:  Bicubic  approximation  (left)and  numerically  integrated  Van  der  Pol  (right)  for  A = 0.1 
with  a 30  by  30  knot  domain 

Table  (3)  shows  the  results  of  changing  the  number  of  knots  keeping  A = 0.1  fixed.  The  third  row 
shows  the  time  T1  which  represents  the  sum  of  the  seconds  it  takes  to  create  the  necessary  matrices 
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in  the  main  program  plus  the  time  it  takes  to  evaluate  once  the  function  called  by  the  optimization 
subroutine.  T2  in  seconds  is  the  sum  of  T1  and  the  time  it  takes  to  run  the  optimization  subroutine 
f solve  in  MATLAB.  The  last  row  gives  the  normed  residual  error  in  each  case. 


Table  2:  Summary  table  for  A = 0.1  showing  timings  and  residual  error. 


A = 0.1 

Knots 

20 

30 

40 

T1  (seconds) 

0 

0.0160 

0.0630 

T2  (seconds) 

7.0620 

13.3290 

38.8750 

Normed  Residual  Error 

1.92718e-010 

9.35666e-010 

6.38028e-009 

Figure  (2)  represents  the  A = 0.6  case.  We  observe  some  deforming  of  the  torus  but  the  bicubic 
surface  appears  to  be  similar  to  that  generated  by  the  MATLAB’s  ode45  solver. 


Figure  2:  Bicubic  approximation  (left)and  numerically  integrated  Van  der  Pol  (right)  for  A = 0.6 
with  a 45  by  45  knot  domain 


Table  3:  Summary  table  for  A = 0.6  showing  timings  and  residual  error. 


A = 0.6 

Knots 

40 

45 

50 

T1  (seconds) 

0.0620 

0.0780 

0.1410 

T2  (seconds) 

68.1250 

80.7190 

265.6710 

Normed  Residual  Error 

8.00207e-011 

9.18135e-010 

7.72686e-008 

As  mentioned  in  the  previous  section,  for  values  of  A > 1.0  the  optimization  did  not  terminate 
successfully  due  to  insufficient  number  of  knots.  Figures  (3)  and  (4)  depict  two  such  cases.  We 
note  that  as  the  strength  of  coupling  increase,  the  inner  sides  of  the  invariant  torus  get  closer  and 
eventually  connect  with  each  other.  Again,  the  results  appear  to  be  consistent. 


7 Conclusions 

This  case  study  presents  an  alternate  approach  to  approximating  the  invariant  torus  of  the  coupled 
Van  der  Pol  oscillator.  The  bicubic  B-spline  collocation  method  in  conjunction  with  the  vectorization 
capabilities  offered  by  the  MATLAB  processing  system  provides  us  with  the  tools  for  developing  a 
significant  time  saving  optimization  algorithm  for  efficiently  computing  the  approximate  periodic 
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Figure  3:  Bicubic  approximation  (left)and  numerically  integrated  Van  der  Pol  (right)  for  A = 1.0 
with  a 50  by  50  knot  domain 


surface.  Comparisons  of  the  present  approach  with  the  MATLAB  numerical  integration  of  the  full 
system  of  coupled  oscillators  indicate  that  the  method  demonstrated  in  this  paper  is  promising.  For 
strong  nolinear  coupling,  computer  memory  appears  to  be  the  only  limitation. 


8 Disclaimer 

Certain  commercial  software  products  are  identified  in  this  paper  in  order  to  adequately  specify  the 
computational  procedures.  Such  identification  does  not  imply  recommendation  or  endorsement  by 
the  National  Institute  of  Standards  and  Technology  nor  does  it  imply  that  the  software  products 
identified  are  necessarily  the  best  available  for  the  purpose. 
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A B-spline  Code  for  Approximate  Torus 

function  Main_VanderPol_Coupled 

°/0MAIN_VAMDERPOL_COUPLED  - This  primary  function  creates  the  diagonal  matrices 

°/0D  and  B for  each  of  the  matrix  equations  approximating  the  functions  r_l  and  r_2 

70and  their  derivatives  at  the  knots  in  both  the  theta  and  phi  directions,  calls  the  subfunctions 

°/0for  computing  the  optimized  B-spline  coefficients  and  plots  the  bicubic 

°/0spline  surfaces. 
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global 

global 

global 

global 

global 

global 

global 

global 

global 

global 

global 


lambda  meul  meu2  a alph 

m n h k deltal  delta2  nm  npl  mpl  nplmpl  mp4  np3  np2  np4  mp3  mp2  np3mp3  np5mp5 

theta  phi  thetavector  phivector 

rl_nm  drldtheta_nm  drldphi_nm 

r2_nm  dr2dtheta_nm  dr2dphi_nm 

thetadot  phidot 

rldot  r2dot 

R1  R2 

Blred  B2red  B3red  Jpat  Jred 
c d 

tim2  tim3 


0/0enter  the  following  parameters 

meul=l . 0 ; °/0input  (’  enter  the  parameter  meul  ,meul=  ’ ) ; 
meu2=l . 414 ; °/0input  (’ enter  the  parameter  meu2  ,meu2=  ’ ) ; 
a=0 . 2 ; %input ( ’ enter  the  parameter  a,a=’); 
alph=0 . 4 ; °/0input  ( ’ enter  the  parameter  alph , alph=  ’ ) ; 


Renter  the  following  inputs 

lambda=input (’ enter  the  parameter  lamda, lambda= ’ ) ; 
n=input (’ enter  harmonics  n,  n=’); 
m=input (’ enter  harmonics  m,  m=’); 


np 1 =n+ 1 ; mp 1 =m+ 1 ; np2=n+2 ; mp2=m+2 ; np3=n+3 ; mp3=m+3 ; 
nplmpl=npl*mpl ; np3mp3=np3*mp3 ; 
nm=n*m ; 


tim=clock ; 

theta=zeros (np3 , 1) ; 
phi=zeros (mp3 , 1 ) ; 

Blred=zeros (nm ,nm) ; 

B2red=zeros(nm,nm) ; 

B3red=zeros (nm , nm) ; 

Jpat=zeros (2*nm , 2*nm) ; 

% mesh  specified  by  {theta=  meul*omegal_i : i=l ,2, . . . ,n+5}&  {phi=meu2*omega2_j : i=l , 2 , . . . ,m+5 

deltal=2*pi/n; 

delta2=2*pi/m ; 


h=2 : np2 ; 

theta (h)=(h-2) *deltal ; 


k=2 :mp2 ; 

phi (k) = (k-2) *delta2 ; 

70creating  the  Blredmatrix 
'/.create  the  matrix  block  B1 
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Bl=zeros (m,m) ; JBl=zeros (m,m) ; 
h=0 ; 

°/0the  first  row  of  the  block 
B1 (rindex(h,0) , r index (h, 0) )=16 ; 

B1 (rindex(h,0) , r index (h, 1) )=4 ; 

B1 (r index (h,0) , r index (h ,m-l) ) =4 ; 

JB1 (rindex(h, 0) ,rindex(h,0) ) =1 ; 

JB1 (rindex(h, 0) ,rindex(h, 1) )=1 ; 

JB1 (rindex(h, 0) ,rindex(h,m-l) ) =1 ; 
7„the  last  row  of  the  block 
B1 (rindex(h,m-l) , r index (h, 0) )=4; 

B1 (rindex(h,m-l) ,r index (h,m-2) )=4; 
B1 (r index (h,m-l) , r index (h,m-l) ) =16 ; 

JB1 (r index (h,m-l) , r index (h , 0) ) = 1 ; 
JB1 (r index (h ,m-l) , r index (h ,m-2) )=1 ; 
JB1 (r index (h,m-l) , r index (h,m-l) )=1 ; 
°/0for  the  intermediate  rows 
for  k=l:m-2; 

B1 (rindex(h,k) ,rindex(h,k) )=16; 

B1 (rindex(h, k) , r index (h,k-l) )=4; 

B1 (r index (h, k) , r index (h ,k+l) ) =4 ; 

JB1 (rindex(h ,k) , r index (h ,k) ) =1 ; 

JB1 (r index (h,k) , r index (h,k-l) )=1 ; 
JB1 (r index (h,k) , r index (h,k+l) )=1 ; 
end 

^create  the  matrix  block  B2 
B2=zeros (m,m) ; JB2=zeros (m,m) ; 
h=0 ; 

°/0the  first  row  of  the  block 
B2 (r index (h,0) , r index (h,0) )=4; 

B2 (r index (h,0) , r index (h, 1) )=1 ; 

B2 (r index (h,0) , r index (h,m-l) )=1 ; 

JB2(rindex(h,0) , rindex(h , 0) ) =1 ; 

JB2 (r index (h,0) , r index (h, 1) )=1 ; 
JB2(rindex(h, 0) , rindex(h ,m-l) )=1 ; 
°/0the  last  row  of  the  block 
B2 (r index (h,m-l) , r index (h,0) )=1 ; 
B2(rindex(h,m-1) ,rindex(h,m-2) )=1 ; 
B2 (r index (h,m-l) , r index (h,m-l) )=4; 

JB2 (r index (h,m-l) , r index (h,0) ) = 1 ; 
JB2(rindex(h,m-l) ,rindex(h,m-2) )=1 ; 
JB2(r index (h,m-l) , r index (h,m-l) )=1 ; 
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for  k=l : m-2 ; 

B2(rindex(h,k) ,rindex(h,k) )=4; 

B2 (r index (h,k) , r index (h , k-1) ) =1 ; 

B2 (r index (h,k) , r index (h,k+l) )=1 ; 

JB2 (rindex(h , k) , r index (h ,k) ) =1 ; 

JB2 (r index (h,k) ,rindex(h,k-l) )=1 ; 

JB2 (rindex(h,k) , rindex(h,k+l) ) = 1 ; 
end 

"/.now  create  the  reduced  matrix  Blred  using  the  above  blocks 
’/.for  the  first  m rows 

Blred(rindex (0 , 0 :m-l) ,rindex(0,0) : rindex(0 ,m-l) ) =B1 ; 

B lred(r index (0 , 0 :m-l) , r index (1 ,0) : rindex(l ,m-l) )=B2 ; 
Blred(rindex(0,0:m-1) ,rindex(m-l ,0) :rindex(m-l ,m-l) )=B2 ; 

Jred (r index (0 , 0 :m-l) , r index (0 , 0) : r index (0 ,m-l) )=JB1 ; 

Jred(r index (0 , 0 : m-1) , r index ( 1 , 0) : r index ( 1 ,m-l) ) =JB2 ; 

Jred (r index (0 , 0 : m-1) , r index (m-1 ,0) : r index (m-1 ,m-l) ) = JB2 ; 

70for  the  last  m rows 

Blred (r index (m-1 , 0 : m-1) , r index (0 , 0) : r index (0 ,m-l) ) =B2 ; 

Blred (r index (m-1 ,0 :m-l) ,r index (m-2 ,0) : r index (m-2 ,m-l) ) =B2 ; 
Blred (r index (m-1 , 0 : m-1) , r index (m-1 , 0) : r index (m-1 , m-1 ) ) =B1 ; 

Jred (r index (m-1 , 0 :m-l) , r index (0 ,0) : r index (0 ,m-l) )=JB2 ; 

Jred (r index (m-1 , 0 : m-1) , r index (m-2 , 0) : r index (m-2 , m-1) ) =JB2 ; 
Jred (r index (m-1 , 0 : m-1) , r index (m-1 , 0) : r index (m-1 , m-1) ) =JB1 ; 

'/.for  the  middle  rows 
for  h=l : m-2 ; 

Blred (r index (h , 0 :m-l) , r index (h , 0 :m-l) )=B1 ; 

B lred(r index (h, 0 : m-1) , r index (h-1 , 0 :m-l) )=B2; 

Blred (r index (h , 0 :m-l) , r index (h+1 , 0 : m-1) ) =B2 ; 

Jred ( r index (h , 0 :m-l) , r index (h, 0 : m-1) ) =JB1 ; 

Jr ed (r index (h,0: m-1) , r index (h-1 , 0 :m-l) ) =JB2 ; 

Jred (r index (h , 0 : m-1) , r index (h+1 , 0 :m-l ) ) =JB2 ; 

end 

clear  B1  B2  JB1  JB2 
%Jacpat 

Jpat (2*nm , 2*nm) =single (0) ; 

Jpat ( 1 : nm , 1 : nm) =Jred ; 

Jpat (1 :nm,nm+l : 2*nm)=Jred; 

Jpat (nm+1 : 2*nm, 1 :nm)=Jred; 

Jpat (nm+1 : 2*nm ,nm+l : 2*nm) =Jred ; 
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clear  Jred; 


"/.creating  the  B2redmatrix 
’/.create  the  matrix  block  B1 
Bl=zeros(m,m) ; 
h=0 ; 

°/0the  first  row  of  the  block 

B1 ( r index (h, 0) , rindex (h, 0) )=12/deltal ; 

B1 (r index (h,0) , r index (h, 1) )=3/deltal ; 

B1 (rindex(h,0) , r index (h,m-l) )=3/deltal ; 

"/.the  last  row  of  the  block 
B1 (rindex(h,m-l) , r index (h, 0) )=3/deltal ; 

B1 (rindex(h,m-l) , r index (h,m-2) )=3/deltal ; 

B1 (rindex(h,m-l) , r index (h,m-l) )=12/deltal ; 

’/.for  the  intermediate  rows 
for  k=l :m-2 ; 

Bl(r index (h,k) , r index (h, k) )=12/deltal ; 

B1 (rindex(h,k) , r index (h , k-1) ) =3/deltal ; 

B1 (rindex(h,k) , r index (h,k+l) )=3/deltal ; 

end 

’/.create  the  matrix  block  B2 
B2=zeros(m,m) ; 
h=0 ; 

°/0the  first  row  of  the  block 

B2 (r index (h , 0) , r index (h,0) )=-12/deltal ; 

B2 (r index (h , 0) , rindex (h, 1) )=-3/deltal ; 

B2(rindex(h,0) , r index (h,m-l) )=-3/deltal ; 

°/0the  last  row  of  the  block 

B2 (r index (h,m-l) , r index (h,0) )=-3/deltal ; 

B2 ( r index (h,m-l) , r index (h,m-2) )=-3/deltal ; 

B2 (rindex(h,m-l) ,rindex(h,m-l) ) =-12/deltal ; 

%for  the  intermediate  rows 
for  k=l:m-2; 

B2 (r index (h,k) , r index (h,k) )=-12/deltal ; 

B2(rindex(h,k) , rindex (h , k-1) )=-3/deltal ; 

B2 (r index (h,k) , r index (h,k+l) )=-3/deltal ; 

end 

'/.now  create  the  reduced  matrix  Bred  using  the  above  blocks 
°/0for  the  first  m rows 

B2red (rindex (0 , 0 :m-l) , rindex ( 1 , 0) : rindex ( 1 ,m-l) ) =B1 ; 
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B2red(rindex(0,0:m-1) ,rindex(m-l ,0) : rindex(m-l ,m-l) )=B2 ; 


%for  the  last  m rows 

B2red(rindex(m-1 ,0 :m-l) ,rindex(0, 0) : rindex(0 ,m-l) )=B1 ; 
B2red(rindex(m-1 ,0 :m-l) ,r index (m-2 , 0) : r index (m-2 ,m-l) )=B2 ; 

°/„for  the  middle  rows 
for  h=l :m-2 ; 

B2red (r index (h,0 :m-l) ,r index (h— 1 ,0:m-l))=B2; 

B2red ( r index (h, 0 :m-l) , r index (h+1 , 0 :m-l) ) =B1 ; 

end 

clear  B1  B2; 

°/0creating  the  B3redmatrix 
%create  the  matrix  block  B1 
Bl=zeros(m,m) ; 
h=0 ; 

°/0the  first  row  of  the  block 
B1 (rindex(h, 0) , r index (h, 1) )=12/delta2; 

B1 (r index (h , 0) , r index (h ,m-l) ) =-12/delta2 ; 

°/othe  last  row  of  the  block 
B1 (rindex(h,m-l) , r index (h,0) )=12/delta2; 

B1 (rindex(h,m-l) , r index (h, m-2) ) =-12/delta2; 

%for  the  intermediate  rows 
for  k=l:m-2; 

B1 (r index (h , k) , r index (h,k-l) ) =-12/ delta2 ; 

B1 (r index (h ,k) , r index (h,k+l) ) =12/delta2 ; 

end 

’/.create  the  matrix  block  B2 

B2=zeros(m,m) ; 

h=0; 

70the  first  row  of  the  block 

B2 (r index (h,0) , r index (h, 1) )=3/delta2; 

B2 (r index (h , 0) , r index (h,m-l) )=-3/delta2; 

%the  last  row  of  the  block 

B2 (rindex(h ,m-l) , r index (h , 0) ) =3/delta2 ; 

B2 ( r index (h,m-l) , r index (h, m-2) )=-3/delta2; 

%for  the  intermediate  rows 
for  k=l :m-2; 

B2 (r index (h , k) , r index (h, k-1 ) ) =-3/delta2 ; 

B2 (r index (h , k) , r index (h,k+l ) ) =3/delta2 ; 

end 
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70now  create  the  reduced  matrix  Bred  using  the  above  blocks 
"/.for  the  first  m rows 
B3red(nm,nm)=single(0) ; 

B3red(r index (0 , 0 :m-l) , r index (0 , 0) : r index (0 ,m-l) ) =B1 ; 

B3red(r index (0 , 0 :m-l) ,r index ( 1 , 0) : r index (1 ,m-l) ) =B2 ; 

B3red(rindex(0 , 0 :m-l) ,rindex(m-l , 0) : rindex(m-l ,m-l) ) =B2 ; 

’/.for  the  last  m rows 

B3red(r index (m-1 ,0 :m-l) ,r index (0,0) : rindex(0 ,m-l) ) =B2 ; 

B3red(rindex(m-1 ,0:m-l) , rindex (m-2 , 0) : rindex(m-2 ,m-l) ) =B2 ; 

B3red(r index (m-1 ,0 :m-l) , rindex (m-1 , 0) : rindex (m-1 ,m-l) )=B1 ; 

‘/.for  the  middle  rows 
for  h=l:m-2; 

B3red (rindex (h , 0 : m-1) , rindex (h, 0 : m-1) )=B1 ; 

B3red (rindex (h , 0 : m-1) , rindex (h-1 , 0 : m-1) ) =B2 ; 

B3red (rindex (h , 0 :m-l) , rindex (h+1 , 0 : m-1) ) =B2 ; 

end 

clear  B1  B2; 

’/.now  for  the  optimization 

“/.initialize  the  parameters  c and  d 
c=zeros (np3mp3 , 1) ; 
d=zeros (np3mp3 , 1 ) ; 

“/.initialize  the  approximation  vectors  rl_nm  & r2_nm 

rl_nm=zeros(nm, 1) ; 

r2_nm=zeros (nm , 1) ; 

phivector=zeros (nm, 1) ; 

thetavector=zeros (nm , 1) ; 

thetadot=zeros (nm, 1) ; phidot=zeros (nm, 1) ; 
rldot=zeros (nm, 1) ; r2dot=zeros (nm, 1) ; 

“/.initialize  the  residuals  R1  and  R2 
Rl=zeros (nm, 1) ; 

R2=zeros (nm, 1) ; 

“/.initialize  the  partial  derivative  vectors  drldtheta  drldphi  dr2dtheta 
“/.dr2dphi 

drldtheta_nm=zeros (nm, 1) ; 
drldphi_nm=zeros (nm, 1) ; 
dr2dtheta_nm=zeros (nm, 1) ; 
dr2dphi_nm=zeros (nm, 1) ; 

“/.initialize  the  matrices  A and  b for  getting  the  initial  conditions 
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bl=zeros  (run,  1)  ; 
b2=zeros(nm, 1) ; 
b=zeros(nm,2) ; 


cred=zeros  (run , 1)  ; 
dred=zeros(nm, 1) ; 

A=zeros  (run , 2)  ; 

A_bar=zeros (run, 2)  ; 

'/(reshape  the  thetal  theta  and  phi  vectors  to  make  each  of  them  a n*m  vector 
phivector=repmat  (phi  (2 : mpl , 1)  ,n , 1)  ; "/.copies  the  whole  vector  n times 
thetavector=kron(theta(2 : npl , 1)  , ones  (m , 1) ) ; "/(copies  each  element  m times 

bl(l :nm,l)=3.53; 
b2 ( 1 : nm , 1) =1 . 18 ; 
cred=inv(Blred) *bl ; 
dred=inv(Blred) *b2 ; 

A= [cred  dred] ; 

°/„load  Jpat 

%’ JacobPattern’ , Jpat , 
timl=etime (clock ,tim) ; 

options=optimset  ( 'Display  ’ , ’ iter  ’ , ’TolX ’ , le-8 , ’LargeScale  ’ , ’ on  ’ , ’ JacobPattern’  , Jpat , ’MaxFunEvals’ 
[A_bar , obj _val] =f solve (Omyobj fun, A, opt  ions) ; 

tim3=etime (clock ,tim) ; 
disp(’Set  up  time  for  F’) 
tim2 

disp( ’Total  time  including  optimization’) 
tim3 

70  use  periodicity  to  get  the  whole  set  of  c and  d 

°/„put  in  the  optimized  c & d values  into  the  cij.dij  vectors 

for  k=0 :n-l ; 

c ( index (k, 0 :m-l) )=A_bar (k*m+l : (k+1) *m, 1) ; 
d ( index (k,0:m-l) )=A_bar (k*m+l : (k+l)*m,2) ; 


end 


°/0put  in  the  periodic  values  of  the  parameters  c 
"/.for  i=0:m-l; 

c ( index (n, 0 :m-l) )=  c (index (0 , 0 : m-1) ) ; 
d ( index (n, 0 :m-l) )=  d ( index (0 , 0 : m-1 ) ) ; 

%end 

’/.for  i=0:n-l; 

c (index (0 :n-l ,m) ) = c (index (0 : n-1 ,0) ) ; 


20 


d ( index (0 : n-1 ,m) ) = d ( index (0 : n-1 , 0) ) ; 
%end 

c ( index (-1 , -1) )=  c( index (n-1 ,m-l) ) ; 

c (index (-1 ,mpl) )=  c( index (n-1 , 1) ) ; 
c (index (npl , -1) )=  c ( index (1 ,m-l) ) ; 

c (index (npl ,mpl) )=  c (index (1,1)); 

d ( index (-1 , -1) )=  d( index (n-1 ,m-l) ) ; 

d ( index (-1 ,mpl) ) = d( index (n-1 , 1) ) ; 
d( index (npl,-l))=  d( index ( 1 ,m-l) ) ; 

d (index (npl ,mpl) )=  d ( index (1 , 1) ) ; 

c (index (n,m) )=  c (index(0 , 0) ) ; 
c (index(0 ,m) )=  c (index(0 , 0) ) ; 

d ( index (n, m) )=  d(index (0 , 0) ) ; 
d(index(0 ,m) )=  d(index (0 , 0) ) ; 

c (index (n-1 ,m) ) = c (index (n-1 , 0) ) ; 
c (index(n ,m-l) ) = c (index(0 ,m-l) ) ; 
c(index(n,mpl))=  c(index(n, 1) ) ; 
c (index (npl ,m) ) = c(index(l ,m) ) ; 

d( index (n-1 ,m) ) = d( index (n-1 , 0) ) ; 
d ( index (n ,m-l) ) = d( index (0 , m-1) ) ; 
d(index(n ,mpl) ) = d(index(n, 1) ) ; 
d(index(npl ,m) ) = d(index(l ,m) ) ; 

'/.for  i=0:m; 

c (index (-1 , 0 :m) ) = c (index (n-1 , 0 :m) ) ; 
c (index (npl , 0 : m) ) = c (index ( 1 , 0 :m) ) ; 

d( index (-1 ,0 :m) )=  d( index (n-1 ,0 :m) ) ; 
d( index (npl ,0 :m) )=  d ( index (1 , 0 :m) ) ; 
‘/.end 

°/0for  i=0:n; 

c (index(0 : n, -1 ) ) = c (index(0 : n,m-l) ) ; 
c (index (0 : n,mpl) ) = c (index (0 : n , 1) ) ; 

d(index(0 :n,-l) )=  d(index(0 : n,m-l) ) ; 
d ( index (0 :n,mpl) )=  d ( index (0 :n, 1) ) ; 
‘/.end 
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°/0plot_surf  ace 

"/(substitute  rapprox  into  thetadot  and  integrate  to  get  thetcap(t) 
tspan=0 : pi/ 100 : 2*pi ; 
thet0=  [0  0] ’ ; 

[t ,thetcap] =ode45(@thetaODEdotBspline_coup,tspan,thetO) ; 
thetlcap=thetcap( : , 1) ; 
thet2cap=thetcap ( : , 2)  ; 

°/0subsitute  thetcap  vector  to  get  the  corresponding  reap  vector 
s=length(thetlcap)  ; 
rlcap=zeros (s , 1) ; 
for  i=l:s; 

rlcap(i)=rlapproxfun_Bspline(thetlcap(i) , thet2cap (i) ) ; 

end 

w=length(thet2cap) ; 
r2cap=zeros (w, 1)  ; 
for  i=l:w; 

r2cap(i)=r2approxfun_Bspline(thetlcap(i) ,thet2cap(i) ) ; 

end 

tspan=0 : pi/100 : 4*pi ; 
z0= [3 . 53  0 1.18  0] ’ ; 

[t  ,zl]  =ode45(@vdpolODEsystemeq_coup,tspan,zO) ; 

rlML=zl  ( : ,1)  ; 

thetalML=zl ( : , 2)  ; 

r2ML=zl ( : ,3) ; 

theta2ML=zl ( : ,4) ; 

/plot  the  full  system  integration  and  the  phase  equation  integration 
numb=input  (’  enter  figure  # numb,  numb^); 
f igure (numb) ; 

title ({’ lambda= ’ .lambda, ,harmonics=’ ,n}) ; 
p=length(thetlcap) ; 
q=length(thet2cap) ; 

x=zeros(p,q) ;y=zeros(p,q) ; z=zeros (p , q) ; 
for  i =l:p 
for  j=l : q 

x(i, j)=r leap (i)  * cos (thetlcap(i) )+r2cap( j ) *cos (thet2cap( j ) ) *cos (thetlcap(i) ) 

y (i , j ) =rlcap(i) *sin(thet lcap(i) )+r2cap( j ) *cos (thet2cap( j ) ) *sin(thetlcap(i) ) 

z(i , j )=r2cap( j ) *sin(thet2cap( j ) ) ; 

end 

end 

subplot (2,2,1) ; 
plot3 (x , y , z)  ; 
grid  on; 
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axis ([-3.5  3.5  -3.5  3.5  -2  2] ) ; 
p=length (thetalML) ; 
q=length(theta2ML) ; 

x=zeros (p , q) ; y=zeros (p , q) ; z=zeros (p , q)  ; 
for  i =l:p 
for  j=l:q 

x(i,j)=rlML(i)*cos (thetalML(i) ) +r2ML( j ) *cos (theta2ML( j ) ) *cos (thetalML (i) ) ; 

y (i , j )=rlML(i) *sin(thetalML(i) )+r2ML(j ) *cos (theta2ML( j ) )*sin (thetalML (i) ) ; 

z(i , j )=r2ML( j ) *sin(theta2ML( j ) ) ; 

end 

end 

subplot (2,2,2) ; 
plot3(x,y ,z) ; 

axis ([-3.5  3.5  -3.5  3.5  -2  2]); 

function  loc=rindex(i , j) 

%RINDEX  - This  function  converts  an  array  matrix  to  a single  column  of  values. 

7. 

70 Input : 

7.  i - matrix  row 

°/0  j - matrix  column 

7. 

"/(Output : 

% loc  - column  vector  of  values 

1 

°/oAUTHOR : 

7.  Sita  Ramamurti 
7.  Mathematics  Program 

Trinity  Washington  University 
7.  125  Michigan  Avenue,  N.E. 

7.  Washington,  D.C.  20017-1094 
7. 

global  m 
loc=i*m+j+l ; 

%************************************************************************** 
l The  Objective  Function 
function  y=myobj f un(A) 

7MY0BJFUN  - This  function  sets  up  the  appropriate  objective  function  to  be 
"/(minimized  using  MATLAB’s  nonlinear  minimization  subroutine. 

7. 

°/o Input : 

7.  A - this  is  the  vector  [c;  d]  of  the  B-spline  coefficients 
7. 

"/oOutput : 
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“/  y - this  is  the  vector  of  "residual  errors" 

% 

“/AUTHOR : 

% Sit a Ramamurti 

7.  Mathematics  Program 
“/  Trinity  Washington  University 
7.  125  Michigan  Avenue,  N.E. 

“/  Washington,  D.C.  20017-1094 

7. 

%************************************************************************** 

global  lambda  meul  meu2  a alph 

global  m n nm  npl  mpl  nplmpl  np3mp3  np5mp5 

global  theta  phi  thetavector  phivector 

global  N M Np  Mp  Npp  Mpp 

global  rl_nm  dr ldtheta_nm  drldphi_nm 

global  r2_nm  dr2dtheta_nm  dr2dphi_nm 

global  rl_lowb  rl_upb  rl_leftb  rl_rightb 

global  r2_lowb  r2_upb  r2_leftb  r2_rightb 

global  rldot_nm  r2dot_nm 

global  thetadot  phidot 

global  rldot  r2dot 

global  R1  R2 

global  Blred  B2red  B3red 
global  tim2 

tim=clock; 

“/  use  the  above  B matrices  to  define  the  approximations  rl_nm  and  r2_nm 
“/  at  the  knots  theta (3)  to  theta (n+3)  and 
% at  phi (3)  to  phi(m+3) 
rl_nm=Blred*A( : , 1) ; 
r2_nm=Blred*A( : ,2) ; 

°/0compute  the  approximations  drldtheta_nm  & dr2dtheta_nm 
drldtheta_nm=B2red* (meul*A( : , 1) ) ; 
dr2dtheta_nm=B2red* (meul*A( : , 2) ) ; 

“/.compute  the  approximations  drldphi_nm  & dr2dphi_nm 
drldphi_nm=B3red* (meu2*A( : , 1) ) ; 
dr2dphi_nm=B3red*(meu2*A( : ,2) ) ; 

“/compute  the  residuals 

thetadot=l+lambda. * (- (1/meul) . * (sin(thetavector) . *cos (thetavector) - . . . 

r l_nm . *s in (thetavector) . "3 . *cos (thetavector) 

a. *r2_nm . *sin (thetavector) . *cos (thetavector) . *sin (phivector) . ~2) ) ; 
phidot =l+lambda . * (- ( l/meu2) . * (sin (phivector) . *cos (phivector) 

alph. *rl_nm. *s in (thetavector) . ~2 . * sin (phivector) . * . . . 
cos (phivector) -r2_nm . * sin (phivector) . “3 . *cos (phivector) ) ) ; 

rldot=lambda. *2 . *rl_nm. * (cos (thetavector) . ~2-rl_nm . * sin (thetavector) . ~2 . *cos (thetavector) . ~2~- 
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a. *r2_nm . *cos (thetavector) . "2 . *sin (phivector) . "2) ; 
r2dot=lambda. *2 . *r2_nm. * (cos (phivector) . ~2-alph. *rl_nm. *sin( thetavector) . ~2 . *cos (phivector) . "2- 
r2_nm. *cos (phivector) . "2 . *sin(phivector) . “2) ; 

"/.compute  the  rl&r2dot  matrices  and  the  residuals 
rldot_nm=drldtheta_nm. *thetadot+drldphi_nm. *phidot ; 
r2dot_nm=dr2dtheta_nm. *thetadot+dr2dphi_nm. *phidot ; 

Rl=rldot_nm-rldot ; 

R2=r2dot_nm-r2dot ; 

y=[Rl;R2] ; 

tim2=etime (clock, tim) ; 
function  loc=index(i , j ) 

“/.INDEX  - This  function  converts  an  array  matrix  to  a single  column  of  values. 

y. 

°/0 Input : 

°/0  i - matrix  row 

"/.  j - matrix  column 

y. 

“/.Output : 

“/.  loc  - column  vector  of  values 

“/. 

“/.AUTHOR: 

“/.  Sit  a Ramamurti 

“/.  Mathematics  Program 

“/.  Trinity  Washington  University 

y.  125  Michigan  Avenue,  N.E. 

“/.  Washington,  D.C.  20017-1094 

% 

global  mp3 
loc=(i+l) *mp3+2+j ; 

“/.integration  of  the  phase  equation  to  get  thetacap 
function  thetaprime=thetaODEdotBspline_coup (t , thet) 

“/.THETAODEDOTBSPLINE  - This  function  substitutes  the  bicubic  spline 

“/.approximations  for  r_l  and  r_2  into  the  Van  der  Pol  system  equations  for  the  derivatives 
°/0of  omega_l  and  omega_2. 

“/.Input : 

“/.  t - a vector  specifying  the  interval  of  integration 

“/.  thet  - a vector  specifying  the  knot  location  in  the  theta  and  phi 

“/.  directions 

“/. 

“/.Output : 

“/.  thetaprime  - vector  of  derivatives  of  omega_l  and  omega_2 
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7. 

"/.AUTHOR: 

"/„  Sita  Ramamurti 

"/„  Mathematics  Program 

"/.  Trinity  Washington  University 

7.  125  Michigan  Avenue,  N.E. 

7.  Washington,  D.C.  20017-1094 

7. 

global  lambda  meul  meu2  a alph 

if  thet(l)<0; 

thet ( 1) =-thet ( 1) ; 

end 

if  thet(2)<0; 

thet (2) =-thet (2) ; 

end 

rl=rlapproxfun_Bspline(thet(l) , thet (2)) ; 
r2=r2approxfun_Bspline(thet(l) , thet (2)) ; 

thetaprime= [1+lambda* (- ( 1/meul) * (sin(meul*thet ( 1) ) *cos (meul*thet (1) ) -rl*sin(meul*thet ( 1) ) '3*. . . 
cos (meul*thet (1) ) - . . . 

a*r2*sin(meul*thet (1) ) *cos (meul*thet (1) ) *sin(meu2*thet (2) ) ~2) ) ; 

1+lambda* (- (1/meul)  * (sin(meu2*thet  (2) ) *cos  (meu2*thet (2) ) -alph*rl*sin(meul*thet(l))"2: 
sin(meu2*thet (2) ) * . . . 

cos (meu2*thet (2) ) -r2*sin(meu2*thet (2) ) "3* cos (meu2*thet (2) ) ) ) ] ; 


‘/.van  der  pol  oscillator  in  polar  coordinates 

%**************  *,*  ********************************************************** 

function  RThetaprime=vdpolODEsystemeq_coup(t ,zl) ; 

°/0VDPOLODESYSTEMEQ_COUP  - This  function  sets  up  the  system  of  coupled  Van 
°/0der  Pol  oscillators  in  polar  form. 

1 

°/0 Input : 

7.  t - a vector  specifying  the  interval  of  integration 
7.  zl  - a vector  of  initial  conditions 

7. 

"/.Output : 

% RTthetaprime  - a column  vector  of  equations  corresponding  to  the  system 

1 

"/.AUTHOR : 

"/„  Sita  Ramamurti 

"/„  Mathematics  Program 

"/.  Trinity  Washington  University 

"/„  125  Michigan  Avenue,  N.E. 
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% Washington,  D.C.  20017-1094 

7. 

global  lambda  meul  meu2  a alph 
RThetaprime= . . . 

[lambda*2*zl (1) * (cos (meul*zl (2) ) ~2-zl (1) *sin(meul*zl (2) ) ~2*cos (meul*zl (2) ) “2- . . . 
a*zl (3) *cos (meul*zl (2) ) ~2*sin(meu2*zl (4) ) “2) ; 

1+lambda* (- (1/meul) ) * (sin (meul *zl (2))*cos (meul*zl (2))-zl(l) * sin (meul *zl (2) ) ~3*cos (meul*zl (2) ) - . . . 

a*zl (3) *sin(meul*zl (2) ) *cos (meul*zl (2) ) *sin(meu2*zl (4) ) ~2) ; 
lambda*2*zl (3) * (cos (meu2*zl (4) ) ~2-alph*zl (1) *sin(meul*zl (2) ) ~2*cos (meu2*zl (4) ) "2- . . . 
zl (3) *cos (meu2*zl (4) ) ~2*sin(meu2*zl (4) ) ~2) ; 

1+lambda* (- (l/meu2) ) * (sin(meu2*zl (4) ) *cos (meu2*zl (4) ) -alph*zl ( 1) *sin(meul*zl (2) ) ~2*sin(meu2*zl (4) ) 
cos (meu2*zl (4) ) -zl (3) *sin(meu2*zl (4) ) ~3*cos (meu2*zl (4) ) ) ] ; 

^********^*****************3f:*3(c3|c^:^:5|e*^:^:**5te********3|c^:*5|c3|c5|c3|c5|e3fc***^:3t:^e3|c5|e3|c5f:5(c*5f:***s|c* 

function  rlhat=r lapproxf un_Bspline (t , s) 

7oRlAPPROXFUN_BSPLINE  - This  function  constructs  the  bicubic  spline 
70surface  r_l  using  the  optimized  spline  coefficients. 

7. 

% Input : 

! t - location  of  the  knot  in  the  theta  direction 

7.  s - location  of  the  knot  in  the  phi  direction 

7. 

XOutput : 

7.  rlhat  - function  value  of  r_l  at  the  (t,s)  knot 

7. 

7.AUTH0R : 

7c  Sita  Ramamurti 
7o  Mathematics  Program 
7c  Trinity  Washington  University 
7c  125  Michigan  Avenue,  N.E. 

7c  Washington,  D.C.  20017-1094 
7c 

global  n deltal  m delta2  np3mp3  np4  mp3  mp4  np3 
global  theta  phi 
global  c 

7oinitialize  the  spline  matrix  B 
Bi=zeros (np3 , 1) ; 

Bj=zeros(mp3, 1) ; 

7»find  out  in  which  interval  ’ t ’ falls 
klo=2;  khi=n+2; 
while  (khi-klo)>l; 

k=f ix( (khi+klo)/2) ; 
if  theta(k)>t; 
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khi=k ; 


else 

klo=k ; 

end 

end 

"/.compute  the  B spline  at  ’t’ 
j =klo ; 
h=deltal ; 

Bi(j-l)=(l/h~3)* (theta ( j+l)-t) '3; 

Bi(j)  =(l/h~3) * (h~3+3*h~2* ( theta ( j + 1) -t) +3*h* (theta ( j + 1) -t) “2-3* (theta( j + 1) -t) ~3) 
Bi ( j+l)=(l/h~3)*(h~3+3*h~2*(t-theta(j) )+3*h* (t-theta( j ) ) ~2-3* (t-theta( j ) ) ~3) ; 
Bi(j+2)=(l/h~3)* (t-theta( j ) ) "3 ; 

"/find  out  in  which  interval  ’s’  falls 
klo=2;  khi=m+2; 
while  (khi-klo)>l; 

k=fix( (khi+klo)/2) ; 
if  phi(k)>s; 
khi=k ; 

else 

klo=k; 

end 

end 

"/compute  the  B spline  at  ’s’ 
v=klo ; 
h=delta2 ; 

Bj(v-l)=(l/h~3)* (phi (v+1) -s) "3 ; 

Bj(v)  =(l/h~3) * (h~3+3*h~2* (phi (v+1) -s) +3*h* (phi (v+1) -s) ~ 2-3* (phi (v+1) -s) ~3) ; 

Bj (v+1) =( l/h~3) * (h~3+3*h~2* (s-phi (v) )+3*h* (s-phi (v) ) ~2-3* (s-phi (v) ) ~3) ; 

Bj (v+2)=(l/h~3) * (s-phi (v) ) "3 ; 

Bivector=zeros (np3mp3 , 1) ; 

Bj vector=zeros (np3mp3 , 1) ; 

Bivector=kron(Bi (1 :np3, 1) , ones (mp3 , 1) ) ; 

B j vector=repmat (Bj ( 1 : mp3 , 1) , np3 , 1) ; 
rlhat=c ’* (Bivector . *Bj vector) ; 

function  r2hat=r2approxf un_Bspline (t ,s) 

°/0R2APPROXFUN_BSPLINE  - This  function  constructs  the  bicubic  spline 
"/surface  r_2  using  the  optimized  spline  coefficients. 

7 

"/Input : 
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"/  t - location  of  the  knot  in  the  theta  direction 

"/  s - location  of  the  knot  in  the  phi  direction 

7. 

°/oOutput : 

"/  r2hat  - function  value  of  r_2  at  the  (t,s)  knot 
"/ 

"/AUTHOR : 

7.  Sita  Ramamurti 
7.  Mathematics  Program 
7.  Trinity  Washington  University 
% 125  Michigan  Avenue,  N.E. 

% Washington,  D.C.  20017-1094 

7. 

global  n deltal  m delta2  np3mp3  np4  mp3  mp4  np3 
global  theta  phi 
global  d 

"/initialize  the  spline  matrix  B 
Bi=zeros (np3 , 1) ; 

Bj=zeros (mp3 , 1) ; 

"/find  out  in  which  interval  ’ t ’ falls 
klo=2;  khi=n+2; 
while  (khi-klo)>l; 

k=f ix ( (khi+klo) /2) ; 
if  theta(k)>t; 
khi=k; 

else 

klo=k; 

end 

end 

"/compute  the  B spline  at  ’ t’ 
j=klo ; 
h=deltal ; 

Bi(j-l)=(l/h~3)*(theta(j+l)-t) ~3; 

Bi(j)  =(l/h~3)* (h~3+3*h~2* (theta ( j +1) -t) +3*h* (theta ( j + 1) -t) ~2-3* (theta( j + 1) -t) ~3) 
Bi(j+1)= (l/h~3) * (h~3+3*h~2* (t-theta( j ) ) +3*h* ( t- theta ( j ) ) “2-3* (t-theta( j ) ) “3) ; 
Bi(j+2)=(l/h~3)* (t-theta( j ) ) ~3 ; 

"/find  out  in  which  interval  ’s’  falls 
klo=2;  khi=m+2; 
while  (khi-klo)>l; 

k=f ix ( (khi+klo) /2) ; 
if  phi(k)>s; 
khi=k ; 
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else 


klo=k; 

end 

end 

/(compute  the  B spline  at  ’s’ 
v=klo ; 
h=delta2 ; 

Bj (v-l)=(l/h~3) * (phi (v+l)-s) ~3; 

Bj(v)  =(l/h~3) * (h'3+3*h'2* (phi (v+1) -s) +3*h* (phi(v+l)-s) '2-3* (phi (v+1) -s) "3) 
Bj(v+l)=(l/h"3)* (h'3+3*h'2* (s-phi (v) )+3*h* (s-phi (v) ) '2-3* (s-phi (v) ) ~3) ; 

Bj (v+2) =(l/h~3)* (s-phi (v)) "3; 

Bivector=zeros (np3mp3 , 1) ; 

B j vector=zeros (np3mp3 , 1 ) ; 

Bivector=kron(Bi (1 : np3 , 1) , ones (mp3 , 1) ) ; 

B j vector=repmat (B j ( 1 : mp3 , 1 ) , np3 , 1 ) ; 
r2hat=d’ * (Bivector . *Bjvector) ; 
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